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DEFLATED  KRYLOV  SUBSPACE  METHODS  FOR  NEARLY 
SINGULAR  LINEAR  SYSTEMS  * 

J.C.  MEZA  AND  W.W.  SYMES  f 

Abstract.  This  paper  concerns  the  use  of  Krylov  subspace  methods  for  the  solution  of  nearly  sin¬ 
gular  nonsymmetric  linear  systems.  We  show  that  the  Incomplete  Orthogonalization  Methods  (IOM) 
in  conjuction  with  certain  deflation  techniques  of  Stewart  and  Chan  can  be  used  to  solve  large  non¬ 
symmetric  linear  systems  which  are  nearly  singular. 

1.  Introduction.  This  study  concerns  the  use  of  Krylov  subspace  methods  for 
the  solution  of  nearly  singular  linear  systems.  There  are  many  problems  in  numerical 
analysis  which  require  the  solution  of  nearly  singular  linear  systems.  For  example, 
in  the  solution  of  nonlinear  systems  by  homotopy  continuation  methods  [1],  or  in 
nonlinear  eigenvalue  problems  [5],  one  often  has  to  solve  nearly  singular  linear  sys¬ 
tems.  Another  example  arises  in  constrained  optimization  problems  [11],  where  the 
constraints  may  be  nearly  linearly  dependent.  Other  examples  include  compartmen- 
tal  models  [10],  and  decomposable  Markov  chains  [16].  We  focus  on  a  problem  from 
seismic  processing  known  as  the  velocity  inversion  'problem  [18]. 

The  goal  of  the  velocity  inversion  problem  is  to  determine  certain  parameters 
describing  the  earth  from  data  taken  at  or  near  the  surface.  The  data  is  known  as 
a  seismogram  and  we  wish  to  determine  the  sound  speed  structure  of  the  earth  from 
these  measurements.  Given  that  we  have  a  functional  relation  defined  by  F{c)  =  s, 
where  c  is  the  speed  of  sound  in  the  earth  and  s  is  a  seismogram,  the  inverse  problem 
is:  given  s  solve  for  c.  Since  the  seismogram  is  usually  contaminated  with  noise  we 
actually  solve  a  nonlinear  least  squares  problem  —  min  J  =  ||s  —  F(c)||  for  a  suitable 
class  of  c. 

A  popular  choice  for  the  solution  of  nonlinear  least  squares  problems  is  the  Gauss- 
Newton  method  [8].  This  method  computes  a  sequence  of  iterates  from  the  formulas 

(1)  (J(ck)T  J(ck))Ac  =  -J(ck)TF(ck), 

(2)  ck+1  =  cfc  +  Ac, 

for  k  =  0, 1, ... ,  starting  from  an  initial  guess  or  model,  c0,  for  the  sound  speed  struc¬ 
ture.  In  the  context  of  the  velocity  inversion  problem  the  elements  of  the  Jacobian 
matrix,  J(cfc),  are  not  readily  available.  However,  it  can  be  shown  that  the  Jacobian 
matrix  acting  on  a  vector  can  be  computed  from  the  solution  of  a  certain  boundary 
value  problem.  Similarly,  the  transpose  of  the  Jacobian  acting  on  a  vector  can  be 
computed  from  the  solution  of  a  related  boundary  value  problem.  The  solutions  of 
the  two  boundary  value  problems  are  independent  which  results  in  a  linear  system 

*  This  work  was  supported  in  part  by  the  National  Science  Foundation  under  grants  DMS-8403148, 
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which  is  almost  symmetric.  The  Gauss-Newton  Hessian,  J(ck)T  J(ck),  would  be  sym¬ 
metric  except  that  the  solution  of  the  two  boundary  value  problems,  which  define 
the  Jacobian  matrix  and  its  transpose  can  only  be  solved  up  to  some  (nonsymmetric) 
discretization  error. 

The  linear  systems  (1)  that  arise  in  the  solution  of  the  velocity  inversion  problem 
by  the  Gauss-Newton  method  are  usually  very  large  and  since  the  matrix  elements 
are  not  readily  available  the  only  recourse  is  to  use  an  iterative  technique  for  the 
solution  of  the  linear  systems.  Additionaly,  due  to  the  physics  of  the  problem  the 
linear  systems  may  be  ill-conditioned  either  by  having  several  large  singular  values 
or  by  having  several  small  singular  values.  The  small  singular  values  arise  because 
the  initial  data  is  band-limited.  Perturbations  in  the  model  problem  corresponding  to 
frequencies  outside  the  passband  of  the  input  data  are  simply  not  seen  by  the  model. 
The  large  singular  values  are  also  inherent  in  this  formulation  of  the  problem.  For 
more  details  the  reader  should  consult  [18]. 

In  this  study  we  address  the  issues  of  computing  the  solution  to  the  linear  sys¬ 
tems  (1)  by  using  an  iterative  technique  which  computes  the  solution  in  the  space 
spanned  by  the  orthogonal  complement  of  the  singular  vectors  corresponding  to  the 
small  singular  values. 

Several  methods  have  been  proposed  for  solving  nearly  singular  linear  systems. 
Consider  the  system  of  linear  equations 

(3)  Ax  —  b, 

where  x  and  b  are  n  dimensional  vectors  and  A  is  an  n  x  n  real  matrix  which  has 
rank  n  —  1.  We  denote  the  set  of  eigenvalues  by  A(A)  =  {A^A),  ...,  A„(A)},  with 
corresponding  eigenvectors  {ui,u2,  The  eigenvalues  are  ordered  so  that 

|Ai|  >  |A2|  >  ...  >  |An|  =  0. 


The  solution  to  (3)  can  be  written  as 


(4) 

ufb 

x  =  xd  +  u 

(5) 

jT 

AW.li 

ii 

B 

where  un  is  a  null  vector  of  A.  The  vector  xd  is  called  the  deflated  solution  to  (3)  and 
(4)  is  called  the  deflated  decomposition.  There  are  many  definitions  of  the  deflated 
solution.  Chan  [3]  defines  deflated  solutions  of  (3)  as  solutions  to  nearby  singular  but 
consistent  systems  derived  from  (3).  The  choice  of  the  nearby  system  will  greatly 
affect  the  deflated  solution.  For  example,  one  might  choose  the  nearest  singular 
matrix  to  A  in  the  Frobenius  norm  and  pick  the  deflated  solution  to  be  the  one  with 
minimum  norm.  It  is  well  known  that  this  choice  amounts  to  setting  the  smallest 
singular  value  of  A  equal  to  zero  in  the  singular  value  decomposition  of  the  matrix  A. 
Other  examples  can  be  found  in  [3]. 
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In  certain  applications  [4]  it  is  preferable  to  compute  the  decomposition  (4).  In 
other  applications,  the  deflated  solution  is  the  only  solution  of  interest.  Notice  that 
if  the  eigenvector  un  were  known  then  both  (4)  and  (5)  could  be  computed  by  first 
computing  x  and  then  orthogonalizing  against  un.  However,  even  if  the  eigenvector 
un  were  known  this  approach  is  not  advisable  because  it  usually  results  in  a  poor 
approximation  to  Xd  due  to  roundoff  errors.  In  particular,  if  the  component  of  the 
solution  in  the  direction  of  the  null  eigenvector  is  large,  then  errors  in  that  component 
tend  to  dominate  the  solution  in  the  other  directions. 

Stewart  [17]  suggested  a  method  for  computing  the  deflated  solution  of  (3)  by  an 
implicit  method.  This  method  uses  orthogonal  projections  constructed  from  approxi¬ 
mations  to  the  singular  vectors  of  the  matrix  A  corresponding  to  the  smallest  singular 
value.  The  disadvantage  of  this  method  is  that  it  requires  a  direct  method  for  the 
solution  of  (3).  Chan  [6]  proposed  a  deflated  Lanczos  method  for  symmetric  positive 
definite  linear  systems  which  only  requires  the  matrix- vector  product  Ax.  The  goal 
of  this  paper  is  to  study  methods  for  computing  the  deflated  solution  for  large  sparse 
nonsymmetric  linear  systems  which  are  nearly  singular. 

The  basic  idea  is  to  use  one  of  the  Krylov  subspace  methods  proposed  by  Saad  [14] 
for  solving  nonsymmetric  linear  systems.  These  methods  compute  an  approximation 
to  the  solution  by  generating  iterates  which  lie  in  a  certain  Krylov  subspace.  The 
approximations  to  the  solution  are  then  computed  by  solving  a  linear  system  described 
by  a  small  upper  Hesssenberg  matrix,  H,  produced  by  Arnoldi’s  method  [2]. 

Arnoldi’s  method  may  be  thought  of  as  a  Galerkin  process  for  approximating 
the  eigenvalues  of  A  by  the  eigenvalues  of  the  upper  Hessenberg  matrix  H.  Like 
the  Lanczos  method,  the  approximations  to  the  eigenvalues  tend  to  be  best  at  the 
extremes  of  the  spectrum,  so  that  the  matrix  H  may  be  expected  to  also  be  nearly 
singular.  Therefore,  one  disadvantage  to  using  the  methods  proposed  by  Saad  is  that 
they  could  require  the  solution  of  a  nearly  singular  linear  system.  Unlike  the  original 
system  (3)  however,  we  need  only  be  able  to  solve  linear  systems  involving  the  much 
smaller  matrix  H .  We  can  then  use  the  deflation  techniques  suggested  by  Stewart 
and  Chan.  In  addition  we  propose  another  method  based  on  solving  a  truncated  least 
squares  problem  analogous  to  the  GMRES  method  suggested  by  Saad  [15]. 

2.  Krylov  Subspace  Methods.  Saad  [14]  proposed  a  class  of  methods  for  solv¬ 
ing  large  sparse  nonsymmetric  linear  systems  based  on  the  Arnoldi  process  [2]  for  com¬ 
puting  the  eigenvalues  of  a  matrix.  Arnoldi’s  method  is  a  generalization  of  the  Lanczos 
method  [13]  for  nonsymmetric  matrices,  and  when  it  is  applied  to  a  symmetric  matrix 
it  reduces  to  the  Lanczos  method.  Like  the  Lanczos  method,  Arnoldi’s  method  is  best 
viewed  as  an  iterative  method  for  approximating  the  eigenvalues  of  large  sparse  matri¬ 
ces.  In  essence,  Arnoldi’s  method  is  just  the  Gram-Schmidt  method  for  computing  an 
orthonormal  basis  for  the  Krylov  subspace  Km(w\,A)  =  span{uq,  Awi, ...,  Am~lW\}. 
The  method  can  by  described  as  follows. 
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Algorithm  2.1.  Arnoldi’s  Method. 

1.  Choose  wi  such  that  ||u>i||  =  1. 

2.  For  j  =  1,2, . . . 

hij  =  ( Awj,wi )  i  =  1,2,...,  j 

i 

wj+ 1  =  Awj  -  ^2  hijWi 

i—i 

Vu  =  IK+ill 

“b'+i  —  Wj+i/hj+ij 


If  we  let  Wm  =  [u>i,  W2, . . . ,  u>m]  ,  then  it  is  easy  to  show  that 
(6)  Hm  -  WlAWm , 

where  the  entries  of  the  upper  Hessenberg  matrix  Hm  are  the  scalars  hij  produced  by 
Arnoldi’s  method  after  m  steps. 

As  Saad  [14]  has  shown,  Arnoldi’s  method  may  be  used  as  a  basis  for  a  class  of 
Krylov  subspace  methods  for  solving  large  sparse  nonsymmetric  linear  systems.  By 
a  Krylov  susbpace  method,  we  mean  any  method  that  approximates  the  solution  to 
the  linear  system  (3)  by  generating  iterates  of  the  form 

xm  —  x0  “H  Znm 

where  x0  is  an  initial  guess,  zm  e  /cm(r0,  A)  and  r0  =  b  —  Ax0.  If  we  carry  out  m 
steps  of  Arnoldi’s  method  starting  with  wx  —  r0/  |[r0||  and  if  we  impose  the  Galerkin 
condition  that  the  residuals  at  each  iteration  be  orthogonal  to  Km(r0,A),  then  this 
yields 


WlAWmym  -  W£r0  =  0. 


Using  the  relation  (6)  yields 


Xm  —  x0  "b  IUotJ/to, 


where  ym  solves  the  system 

(7)  Hmym  =  ||r0|l  eu 

and  ei  =  (1,  0, ...,  0)T. 

This  defines  the  Full  Orthogonalization  Method  [14],  for  solving  (3). 
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ALGORITHM  2.2.  Full  Orthogonalization  Method. 

1.  Choose  £0  and  compute  r0  =  h  —  Ax0.  Set  Wi  =  r0/||r0||. 

2.  For  j  =  1,2 ,  ...m 

h{j 
™j+ 1 

™j+i 

3.  Form  the  solution: 

•  Solve  tfmym  =  Urolle! 

•  Set  =  a:0  +  Wmym. 

In  practice,  the  number  of  iterations  is  chosen  so  that  the  approximate  solution  xm 
is  sufficiently  accurate.  Usually  this  is  measured  by  requiring  that  the  initial  residual 
be  decreased  by  a  user-specified  amount.  Fortunately,  the  residual  at  any  iteration 
may  be  computed  without  actually  computing  the  solution  to  (3)  through  the  relation 

[14], 

(®)  11^  •/4'cm||  =  ^m+l,7n|emym|- 

Although  the  computation  of  the  residual  by  (8)  requires  solving  (7)  for  ym  there 
are  ways  to  circumvent  this  computation  by  carrying  an  LU  or  QR  factorization  of 
H  throughout  the  Arnoldi  process.  If  after  m  iterations  the  approximate  solution  has 
not  converged  then  it  is  possible  to  restart  the  algorithm  using  the  current  estimate 
of  x  as  the  new  initial  guess.  This  method  is  denoted  by  FOM(fc),  or  the  restarted 
FOM. 

It  is  well  known  that  the  Arnoldi  process  may  be  viewed  as  a  Galerkin  process  for 
estimating  the  eigenvalues  of  a  matrix.  In  particular,  if  we  apply  Algorithm  2.2  to  a 
linear  system  that  is  nearly  singular  then  the  upper  Hessenberg  matrix,  Hm,  which  is 
generated  after  m  steps  of  the  Arnoldi  process  will  probably  have  a  small  eigenvalue. 
Therefore  if  we  solve  (7)  for  ym  in  the  straightforward  way,  our  computed  solution 
will  be  inaccurate  for  the  reasons  indicated  in  Section  1 . 

Fortunately,  computing  the  deflated  solution  of  (7)  is  easier  than  computing  the 
deflated  solution  of  (3).  Since  the  matrix  Hm  has  dimension  m  <C  n  the  solution  of 
(7)  is  at  least  computationally  easier.  Moreover,  the  matrix  elements  of  Hm  are  on 
hand  whereas  the  matrix  elements  of  A  are  not  available.  The  next  section  describes 
some  of  the  deflation  techniques  which  can  be  used  to  compute  the  deflated  solution 
of  (7)  in  a  stable  manner. 

3.  Deflation  Methods.  This  section  describes  three  methods  for  computing 
the  deflated  solution  to  a  nearly  singular  linear  system.  In  particular,  we  consider  the 


=  (Aw^w,)  i  —  1, 2,  ...,j 
i 

-  Awj  -  ^2  hijWi 
i=i 

=  \\™j+i\\ 
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linear  system 

(9) 


Hmy  =  /, 


where  y  and  /  are  m  dimensional  vectors  and  Hm  is  the  upper  Hessenberg  matrix 
generated  after  m  steps  of  Arnoldi’s  method.  Further,  we  assume  that  Hm  is  nearly 
singular  with  a  rank  deficiency  of  at  most  one.  The  extensions  to  null  spaces  of 
dimensions  greater  than  one  will  be  treated  in  a  later  section. 

The  first  method,  which  uses  orthogonal  projections,  was  proposed  by  Stewart 
[17].  The  second  method  is  a  generalization  of  a  technique  suggested  by  Chan  [6]  for 
symmetric  positive  definite  systems.  In  essence,  this  method  uses  a  QR  iteration  to 
decouple  the  linear  system  into  a  well  conditioned  problem  plus  a  component  corre¬ 
sponding  to  the  small  eigenvalue.  The  third  method  computes  the  deflated  solution 
by  solving  a  truncated  least  squares  problem. 

3.1.  Deflation  by  Orthogonal  Projection.  Stewart  [17]  proposed  an  implicit 
method  for  computing  the  deflated  solution  to  a  nearly  singular  linear  system.  His  al¬ 
gorithm  consists  of  constructing  two  orthogonal  projectors  defined  by  approximations 
to  the  singular  vectors  corresponding  to  the  small  singular  value. 

Consider  the  right  and  left  singular  vectors,  respectively  v  and  u,  corresponding 
to  the  smallest  singular  value  of  Hm.  Define  the  orthogonal  projectors  Pu  —  I  —  uuT 
and  Pv  —  I  —  vvT .  The  projector  PU(PV)  is  merely  the  orthogonal  projector  onto  the 
orthogonal  complement  of  the  space  spanned  by  u(v).  It  is  easy  to  show  that  the 
deflated  solution  to  (9)  is  the  unique  vector  satisfying  the  relations 

(10)  PuHmPvyd  =  Puf , 

(n)  PvVd  =  yd- 

Stewart  suggests  using  the  following  algorithm  based  on  iterative  refinement  to  solve 
for  yd  ■ 

ALGORITHM  3.1.  Deflation  by  Orthogonal  Projection. 

1.  Set  y  =  0 

2.  For  k  =  1,2, ... 

•  Solve  Hmd  =  Pu(f  -  Hmy ), 

•  Set  y  =  y  +  Pvd. 

3.  yd  =  y 

Since  the  singular  vectors  are  not  known,  Stewart  suggests  approximating  u  and  v 
by  a  variant  of  the  inverse  power  method.  In  the  case  of  interest,  where  the  small  sin¬ 
gular  value  is  isolated  the  inverse  power  method  is  known  to  converge  rapidly.  Stewart 
also  gives  conditions,  which  depend  on  the  accuracy  attained  in  the  approximation  to 
the  vectors  u  and  v,  under  which  Algorithm  3.1  will  converge  to  the  deflated  solution 
of  (9). 
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3.2.  Deflation  by  QR  iteration.  The  second  method  for  computing  the  de¬ 
flated  solution  is  a  generalization  of  a  method  proposed  by  Chan  [6]  for  symmetric 
positive  definite  systems.  To  motivate  the  discussion,  assume  that  we  have  the  pair 
( Am , um  )  of  the  unreduced  upper  Hessenberg  matrix  Hm.  It  is  well  known  that  one 
step  of  the  shifted  QR  iteration  method  with  a  shift  of  \m  reduces  the  matrix  Hm  so 
that  the  eigenvalue  Am  appears  on  the  diagonal  and  the  corresponding  subdiagonal 
element  is  zero.  The  linear  system  (9)  can  then  be  decoupled  so  that  it  is  easily  solved 
for  the  deflated  solution. 

In  particular,  if  we  compute  the  matrix  H M  by  the  following  two  step  procedure: 
Compute  the  QR  factorization  of Hm  -  Ara/Form  H W  =  RQ  +  A mJ,  then  the  matrix 
H ^  takes  on  the  form 


m  —  1  1 

m  —  1  /  H  h 
1  V  0  Am 


where  H  is  an  upper  Hessenberg  matrix  of  order  m  —  1.  It  is  easy  to  show  that 
#(1)  =  QTHmQ  ,  so  that  we  can  transform  (9)  into 


(12) 


=  /, 


where  y  =  Qz  and  /  =  QTf. 

If  we  partition  the  vectors  z  and  /  so  that 


z  = 


22 


then  the  deflated  solution  to  (12)  is 

(13)  zd  = 


/ 


E-'h 

0 


fx 

h 


The  deflated  solution  to  (9)  can  now  be  computed  from  yd  —  Qzd. 

There  are  two  remarks  in  order.  The  first  is  that  in  general  we  do  not  know  the 
eigenvalue  Am;  however,  as  in  the  previous  section,  we  may  compute  an  approximation 
to  Xm  by  using  the  inverse  power  method.  The  second  remark  is  that  even  if  we  had 
an  exact  value  for,  Am,  in  practice  the  matrix  Hm  will  not  be  reduced  in  one  QR 
iteration  step.  Wilkinson  [19]  suggests  iterating  with  the  shifted  QR  method  until  the 
element  on  the  subdiagonal  has  converged  to  zero.  We  use  the  test 

(W)  <  e 

to  check  for  convergence.  Two  or  three  iterations  are  usually  sufficient  to  satisfy  (14). 
This  leads  to  the  following  algorithm. 
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Algorithm  3.2.  Deflation  by  QR  iteration. 

1.  Compute  Xm  from  Hm  via  inverse  iteration. 

2.  tf(°)  =  Hm 

3.  For  j  =0,1,...  until  convergence. 

•  Compute  the  QR  factorization  of  —  \mI. 

•  Compute  =  RQ  +  A mI. 

4.  Set  zd  =  [^-1/1,0]T. 

5.  Set  yd  =  Qzd. 

If  the  dimension  of  the  null  space  is  greater  than  one  then  the  issues  become 
more  complicated.  There  are  two  cases  to  consider:  a  small  eigenvalue  of  multiplicity 
greater  than  one  and  the  case  of  several  distinct  small  eigenvalues.  The  case  of  a  small 
real  eigenvalue  with  a  multiplicity  greater  than  one  is  easily  handled  since  the  shifted 
QR  method  will  still  converge.  The  case  of  several  distinct  small  eigenvalues  is  harder 
to  address  since  the  QR  method  does  not  guarantee  that  the  converged  eigenvalues 
will  appear  in  any  particular  order  on  the  diagonal  of  A  related  issue  is  that 

of  complex  eigenvalues.  Since  the  inverse  power  method  converges  to  the  eigenvalue 
of  smallest  modulus,  we  must  be  careful  in  choosing  the  shift.  A  good  choice  would 
be  to  use  Francis’  implicit  double  shift  QR  method  and  with  a  shift  as  described  by 
Wilkinson  [19]. 

3.3.  Deflation  by  Truncated  Least  Squares.  The  last  method  we  discuss  is 
derived  from  a  technique  for  solving  rank  deficient  least  squares  problems  [12,  pp. 
162-167].  The  general  idea  in  these  problems  is  to  compute  a  QR  factorization  of  the 
nearly  singular  matrix  Hm  such  that  the  elements  on  the  diagonal  of  the  matrix  R 
display  the  rank  deficiency.  The  solution  to  the  linear  system  (9)  is  then  computed 
by  setting  the  small  diagonal  elements  of  the  matrix  R  equal  to  zero  and  solving  a 
truncated  least  squares  problem. 

Assume  that  the  matrix  Hm  has  exactly  one  zero  eigenvalue  and  consider  its  QR 
factorization, 

HmIl  =  QR, 

where  II  is  a  permutation  matrix  chosen  so  that  the  matrix  R  has  the  form 

m  —  1  1 

_  m  —  if  Ru  R12  A 

R=  1  [  0  0  )' 

where  Rn  is  upper  triangular. 

In  this  case,  the  least  squares  solution  to  (9)  can  be  easily  computed.  In  fact, 

\\Hmy  —  f\\2  —  H-Rn^i  —  (A  —  ^12-^2)  ||2  +  llAll2? 

where 

(15)  n  Ty=\Zl 
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and 


(is)  qti=  i . 

If  we  set  z2  —  0,  then  the  solution  is  given  by 

(u)  «  =  n  Rn0fl  . 

The  vector  y  is  called  the  basic  solution.  In  general,  it  is  not  the  least  norm  least 
squares  solution  unless  R12  =  0. 

In  practice,  the  matrix  Hm  will  never  have  an  eigenvalue  exactly  equal  to  zero. 
The  question  is  then  to  determine  the  rank  of  the  matrix  Hm  from  the  elements  of  R. 
One  popular  choice  is  the  method  of  column  pivoting  implemented  in  LINPACK  [9]. 
This  method  is  usually  reliable  in  detecting  the  rank  deficiency  of  a  matrix,  although 
there  are  counterexamples  where  the  method  may  fail.  Chan  [7]  suggested  another 
method  for  producing  a  QR  factorization  which  guarantees  a  small  R22  element.  The 
essential  ideas  can  be  explained  in  the  case  of  a  rank  one  deficient  matrix.  First,  we 
need  the  following  theorem  proved  by  Chan  [7]. 

Theorem  1.  Suppose  that  x  6  IRn  with  ||x||  =  1  such  that  ||Aa;||  =  e.  Let  H  be  a 
permutation  such  that  if  IITa;  =  y,  then  \yn\  =  HyH^  .  Then  if  AH  =  QR  is  the  QR 
factorization  of  AH,  then 

Knl  <y/ne. 


The  usefulness  of  this  theorem  is  apparent  if  we  consider  the  right  singular  vector, 
v,  of  the  matrix  Hm  corresponding  to  the  smallest  singular  value  am.  Then  we  have 
H  =  1,  and  \\Hmv\\  =  am.  If  we  define  the  permutation  II  by 

=  IMU 

then  HmH  =  QR  has  a  pivot  rmm  at  least  as  small  as  y/mam  in  absolute  value.  As 
in  Stewart’s  method,  all  that  is  required  is  an  approximation  to  v,  which  may  be 
computed  by  the  inverse  power  method.  This  suggests  the  following  algorithm  for 
solving  (9). 

ALGORITHM  3.3.  Deflation  by  Truncated,  Least  Squares. 

1.  Compute  the  QR  factorization  of  Hm. 

2.  Compute  v  and  am  from  Hm  via  inverse  power  method. 

3.  Compute  II  so  that  |(IIri;)m|  =  HuH^ 

4.  Compute  the  QR  factorization  of  HmH. 

5.  Compute  y  =  II  i?n/i,0]T. 


9 


This  method  has  the  obvious  advantage  of  being  immediately  applicable  to  null 
spaces  with  a  dimension  greater  than  one.  However  the  question  of  rank  deficiency  is 
still  a  hard  problem  and  the  user  must  be  able  to  supply  a  tolerance  which  specifies 
the  amount  of  ill- conditioning  allowed. 


4.  Numerical  Results..  This  section  presents  several  numerical  experiments 
comparing  the  various  methods  described  in  Section  3. 

Recall  that  he  linear  system  of  interest  is 


Ax  =  b , 


where  x  and  b  are  n  dimensional  vectors  and  A  is  an  n  x  n  real  matrix  which  is  nearly 
singular.  All  of  the  numerical  results  presented  here  are  for  linear  systems  of  order 
10.  The  results  are  similar  for  larger  systems.  The  numerical  experiments  were  run 
on  a  Sun  3/160  computer,  using  single  precision  arithmetic  (machine  epsilon  «  10-7). 
The  method  was  said  to  converge  whenever 


M 

INI 


<  10 


-6 


The  basic  method  used  was  the  restarted  FOM(k).  To  avoid  the  issue  of  deciding 
when  the  eigenvalues  of  the  upper  Hessenberg  matrix,  H ,  are  good  approximations 
to  the  eigenvalues  of  the  matrix  A,  we  set  k  =  n. 

The  first  test  case  was  a  small  peturbation  to  a  symmetric  positive  definite  matrix. 
Define 


A(e)  =  D  +  eE. 

The  matrix  D  is  defined  by  D  =  diag  (10_7,2,3, .. .  ,  n),  and  I  varies  from  1  to  7  . 
The  nonsymmetric  perturbations  are  generated  using  the  random  number  generator, 
URAND,  from  IMSL.  The  matrices,  E ,  are  computed  by  generating  uniform  random 
numbers  between  [-0.5,  +0.5  ],  and  normalizing  so  that  ||£?||2  =  1.  The  amount  of 
nonsymmetry  can  then  be  adjusted  by  varying  e.  We  set  e  =  10~3. 

The  second  set  of  test  cases  was  picked  from  a  study  done  by  Chan  [6].  In  this 
set  of  problems  we  set 

A  =  (J  —  uuT)  D  (/  —  vvT), 

where  u  and  v  were  chosen  randomly  with  the  constraint  that  they  have  norm  one. 
The  matrix  D  was  chosen  as  in  problem  1. 

The  purpose  of  the  third  problem  is  to  simulate  a  typical  linear  system  arising  in 
the  velocity  inversion  problem.  These  problems  usually  have  one  or  more  small  sin¬ 
gular  values  and  one  or  more  large  singular  values  with  the  rest  of  the  spectrum  fairly 
well-conditioned.  In  this  problem,  we  set  the  matrix  D  =  diag  (10_/,  1, . . . ,  3, 3000)  , 
with  the  values  of  d2  through  dn_i  varying  uniformly  between  1  and  3.  The  matrix  A 
is  then  computed  as  in  problem  1.  This  example  generates  a  well- conditioned  problem 
if  the  small  and  large  eigenvalues  are  excluded,  which  is  typical  of  some  of  the  velocity 
inversion  problems. 

The  four  methods  discussed  in  Section  3  are: 
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1.  Full  Orthogonalization  Method  with  deflation  after 
the  solution  of  Ax  =  b.  —  FOM. 

2.  FOM  with  Stewart’s  orthogonal  projection  deflation  —  FOMOP. 

3.  FOM  with  the  truncated  least  squares  deflation  —  FOMLSQ. 

4.  FOM  with  the  QR  iteration  deflation  —  FOMQRI. 

Figures  1-3  plot  the  log||r||  versus  I.  Since  we  are  interested  in  the  deflated 
solution  we  have  chosen  to  plot  the  norm  of  the  residual  corresponding  to  the  deflated 
solution.  As  the  value  of  I  increases,  the  linear  systems  become  more  ill-conditioned, 
with  the  smallest  singular  value  of  each  matrix  approximately  equal  to  10_/.  The  effect 
of  the  near  singularity  of  the  linear  systems  on  Method  1  is  apparent  as  the  plots  show 
the  norm  of  the  residual  increasing  as  the  linear  system  becomes  more  ill-conditioned. 
This  is  to  be  expected  as  the  error  in  the  component  corresponding  to  the  smallest, 
singular  values  tends  to  contaminate  the  rest  of  the  solution.  Method  2,  FOMOP, 
is  clearly  the  best  method  in  terms  of  producing  the  smallest  residuals.  Method  3, 
FOMLSQ,  was  disappointing  in  that  the  norm  of  the  residual  was  consistently  worse 
than  the  other  methods.  Method  4,  FOMQRI,  was  inconsistent  in  this  set  of  problems. 
In  problems  1  and  3,  FOMQRI  was  almost  as  good  as  FOMOP.  However  in  problem 
2  the  method  performed  much  worse. 

The  question  of  extending  these  methods  to  null  spaces  of  higher  dimensions  is 
also  of  interest.  In  this  respect,  FOMLSQ  can  be  easily  extended  while  the  other 
methods  would  require  some  extra  work.  In  fact  the  ease  with  which  FOMLSQ  can 
be  extended  to  solving  systems  which  have  several  small  singular  values  is  perhaps 
the  only  redeeming  factor  of  this  method.  Figures  4-6  illustrate  the  effect  of  solving 
the  same  systems  of  problems  1-3  using  both  a  truncated  least  squares  approach  and 
simply  solving  the  system  using  a  full  least  square  problem.  It  can  be  seen  that  the  two 
curves  usually  intersect  between  1  —  2  and  7  =  3  which  corresponds  to  k(A)  fa  1000. 
This  implies  that  using  FOMLSQ  with  a  properly  chosen  tolerance  would  perform 
better  than  always  solving  the  truncated  least  squares  problem.  Fortunately  this  is 
the  case  of  interest.  This  approach  would  still  not  be  as  good  as  using  FOMOP,  but 
it  has  the  advantage  that  it  is  easier  to  extend  to  null  spaces  with  dimension  greater 
than  one,  whereas  extending  FOMOP  to  handle  several  small  singular  values  would 
require  extra  work. 

Further  experiments  remain  to  be  done.  It  is  not  clear  which  deflated  solution 
is  best  in  terms  of  the  norm  of  the  error.  In  fact,  there  are  many  definitions  of  the 
deflated  solution  and  the  choice  of  the  deflated  solution  greatly  alters  the  results. 
That  is  one  reason  we  have  chosen  to  work  with  the  norm  of  the  residual.  Another 
possibility  lies  in  solving  the  upper  Hessenberg  system  by  using  the  singular  value 
decomposition.  Although  this  may  seem  like  to  much  work  at  first  sight,  the  prob¬ 
lems  we  are  dealing  with  have  the  property  that  one  matrix-vector  multiply  is  very 
expensive.  In  this  context,  computing  the  SVD  for  a  small  upper  Hessenberg  matrix 
would  be  insignificant.  This  approach  also  has  the  added  attraction  that  (like  the 
FOMLSQ  method)  it  is  easily  extended  to  null  spaces  of  dimension  greater  than  one. 
This  research  will  be  the  subject  of  a  later  report. 
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Figure  1.  Problem  1  -  Log 


=FOM,  2=F0M0P,  3=F0MLSQ,  4=F0MQRI 


Figure  3.  Problem  3  -  Log 


FOM,  2=F0M0P,  3=F0MLSQ,  4=F0MQRI 
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